home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / fft / signals_source.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-05-14  |  6.4 KB  |  289 lines
  1. /* fft/signals_source.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include "signals.h"
  21.  
  22. int
  23. FUNCTION(fft_signal,complex_pulse) (const size_t k,
  24.                     const size_t n,
  25.                     const size_t stride,
  26.                     const BASE z_real,
  27.                     const BASE z_imag,
  28.                     BASE data[],
  29.                     BASE fft[])
  30. {
  31.   size_t j;
  32.  
  33.   if (n == 0)
  34.     {
  35.       GSL_ERROR ("length n must be positive integer", GSL_EDOM);
  36.     }
  37.  
  38.   /* gsl_complex pulse at position k,  data[j] = z * delta_{jk} */
  39.  
  40.   for (j = 0; j < n; j++)
  41.     {
  42.       REAL(data,stride,j) = 0.0;
  43.       IMAG(data,stride,j) = 0.0;
  44.     }
  45.  
  46.   REAL(data,stride,k % n) = z_real;
  47.   IMAG(data,stride,k % n) = z_imag;
  48.  
  49.   /* fourier transform, fft[j] = z * exp(-2 pi i j k / n) */
  50.  
  51.   for (j = 0; j < n; j++)
  52.     {
  53.       const double arg = -2 * M_PI * ((double) ((j * k) % n)) / ((double) n);
  54.       const BASE w_real = (BASE)cos (arg);
  55.       const BASE w_imag = (BASE)sin (arg);
  56.       REAL(fft,stride,j) = w_real * z_real - w_imag * z_imag;
  57.       IMAG(fft,stride,j) = w_real * z_imag + w_imag * z_real;
  58.     }
  59.  
  60.   return 0;
  61.  
  62. }
  63.  
  64.  
  65. int
  66. FUNCTION(fft_signal,complex_constant) (const size_t n,
  67.                        const size_t stride,
  68.                        const BASE z_real,
  69.                        const BASE z_imag,
  70.                        BASE data[],
  71.                        BASE fft[])
  72. {
  73.   size_t j;
  74.  
  75.   if (n == 0)
  76.     {
  77.       GSL_ERROR ("length n must be positive integer", GSL_EDOM);
  78.     }
  79.  
  80.   /* constant, data[j] = z */
  81.  
  82.   for (j = 0; j < n; j++)
  83.     {
  84.       REAL(data,stride,j) = z_real;
  85.       IMAG(data,stride,j) = z_imag;
  86.     }
  87.  
  88.   /* fourier transform, fft[j] = n z delta_{j0} */
  89.  
  90.   for (j = 0; j < n; j++)
  91.     {
  92.       REAL(fft,stride,j) = 0.0;
  93.       IMAG(fft,stride,j) = 0.0;
  94.     }
  95.  
  96.   REAL(fft,stride,0) = ((BASE) n) * z_real;
  97.   IMAG(fft,stride,0) = ((BASE) n) * z_imag;
  98.  
  99.   return 0;
  100.  
  101. }
  102.  
  103.  
  104. int
  105. FUNCTION(fft_signal,complex_exp) (const int k,
  106.                   const size_t n,
  107.                   const size_t stride,
  108.                   const BASE z_real,
  109.                   const BASE z_imag,
  110.                   BASE data[],
  111.                   BASE fft[])
  112. {
  113.   size_t j;
  114.  
  115.   if (n == 0)
  116.     {
  117.       GSL_ERROR ("length n must be positive integer", GSL_EDOM);
  118.     }
  119.  
  120.   /* exponential,  data[j] = z * exp(2 pi i j k) */
  121.  
  122.   for (j = 0; j < n; j++)
  123.     {
  124.       const double arg = 2 * M_PI * ((double) ((j * k) % n)) / ((double) n);
  125.       const BASE w_real = (BASE)cos (arg);
  126.       const BASE w_imag = (BASE)sin (arg);
  127.       REAL(data,stride,j) = w_real * z_real - w_imag * z_imag;
  128.       IMAG(data,stride,j) = w_real * z_imag + w_imag * z_real;
  129.     }
  130.  
  131.   /* fourier transform, fft[j] = z * delta{(j - k),0} */
  132.  
  133.   for (j = 0; j < n; j++)
  134.     {
  135.       REAL(fft,stride,j) = 0.0;
  136.       IMAG(fft,stride,j) = 0.0;
  137.     }
  138.  
  139.   {
  140.     int freq;
  141.  
  142.     if (k <= 0)
  143.       {
  144.     freq = (n-k) % n ;
  145.       }
  146.     else
  147.       {
  148.     freq = (k % n);
  149.       };
  150.  
  151.     REAL(fft,stride,freq) = ((BASE) n) * z_real;
  152.     IMAG(fft,stride,freq) = ((BASE) n) * z_imag;
  153.   }
  154.  
  155.   return 0;
  156.  
  157. }
  158.  
  159.  
  160. int
  161. FUNCTION(fft_signal,complex_exppair) (const int k1,
  162.                       const int k2,
  163.                       const size_t n,
  164.                       const size_t stride,
  165.                       const BASE z1_real,
  166.                       const BASE z1_imag,
  167.                       const BASE z2_real,
  168.                       const BASE z2_imag,
  169.                       BASE data[],
  170.                       BASE fft[])
  171. {
  172.   size_t j;
  173.  
  174.   if (n == 0)
  175.     {
  176.       GSL_ERROR ("length n must be positive integer", GSL_EDOM);
  177.     }
  178.  
  179.   /* exponential,  data[j] = z1 * exp(2 pi i j k1) + z2 * exp(2 pi i j k2) */
  180.  
  181.   for (j = 0; j < n; j++)
  182.     {
  183.       const double arg1 = 2 * M_PI * ((double) ((j * k1) % n)) / ((double) n);
  184.       const BASE w1_real = (BASE)cos (arg1);
  185.       const BASE w1_imag = (BASE)sin (arg1);
  186.       const double arg2 = 2 * M_PI * ((double) ((j * k2) % n)) / ((double) n);
  187.       const BASE w2_real = (BASE)cos (arg2);
  188.       const BASE w2_imag = (BASE)sin (arg2);
  189.       REAL(data,stride,j) = w1_real * z1_real - w1_imag * z1_imag;
  190.       IMAG(data,stride,j) = w1_real * z1_imag + w1_imag * z1_real;
  191.       REAL(data,stride,j) += w2_real * z2_real - w2_imag * z2_imag;
  192.       IMAG(data,stride,j) += w2_real * z2_imag + w2_imag * z2_real;
  193.     }
  194.  
  195.   /* fourier transform, fft[j] = z1 * delta{(j - k1),0} + z2 *
  196.      delta{(j - k2,0)} */
  197.  
  198.   for (j = 0; j < n; j++)
  199.     {
  200.       REAL(fft,stride,j) = 0.0;
  201.       IMAG(fft,stride,j) = 0.0;
  202.     }
  203.  
  204.   {
  205.     int freq1, freq2;
  206.  
  207.     if (k1 <= 0)
  208.       {
  209.     freq1 = (n - k1) % n;
  210.       }
  211.     else
  212.       {
  213.     freq1 = (k1 % n);
  214.       };
  215.  
  216.     if (k2 <= 0)
  217.       {
  218.     freq2 = (n - k2) % n;
  219.       }
  220.     else
  221.       {
  222.     freq2 = (k2 % n);
  223.       };
  224.  
  225.     REAL(fft,stride,freq1) += ((BASE) n) * z1_real;
  226.     IMAG(fft,stride,freq1) += ((BASE) n) * z1_imag;
  227.     REAL(fft,stride,freq2) += ((BASE) n) * z2_real;
  228.     IMAG(fft,stride,freq2) += ((BASE) n) * z2_imag;
  229.   }
  230.  
  231.   return 0;
  232.  
  233. }
  234.  
  235.  
  236. int
  237. FUNCTION(fft_signal,complex_noise) (const size_t n,
  238.                     const size_t stride,
  239.                     BASE data[],
  240.                     BASE fft[])
  241. {
  242.   size_t i;
  243.   int status;
  244.  
  245.   if (n == 0)
  246.     {
  247.       GSL_ERROR ("length n must be positive integer", GSL_EDOM);
  248.     }
  249.  
  250.   for (i = 0; i < n; i++)
  251.     {
  252.       REAL(data,stride,i) = (BASE)urand();
  253.       IMAG(data,stride,i) = (BASE)urand();
  254.     }
  255.  
  256.   /* compute the dft */
  257.   status = FUNCTION(gsl_dft_complex,forward) (data, stride, n, fft);
  258.  
  259.   return status;
  260. }
  261.  
  262.  
  263. int
  264. FUNCTION(fft_signal,real_noise) (const size_t n,
  265.                  const size_t stride,
  266.                  BASE data[],
  267.                  BASE fft[])
  268. {
  269.   size_t i;
  270.   int status;
  271.  
  272.   if (n == 0)
  273.     {
  274.       GSL_ERROR ("length n must be positive integer", GSL_EDOM);
  275.     }
  276.  
  277.   for (i = 0; i < n; i++)
  278.     {
  279.       REAL(data,stride,i) = (BASE)urand();
  280.       IMAG(data,stride,i) = 0.0;
  281.     }
  282.  
  283.   /* compute the dft */
  284.   status = FUNCTION(gsl_dft_complex,forward) (data, stride, n, fft);
  285.  
  286.   return status;
  287. }
  288.  
  289.